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ABSTRACT 

We revisit the issue of sensitivity to initial flow and intrinsic variability in hot-Jupiter atmo¬ 
spheric flow simulations, originally investigated by Cho et al[| ( |2008[ ) and Thrastarson & Cho 
(20101. The flow in the lower region (~1 to 20 MPa) ‘dragged’ to immobility and uniform 


temperature on a very short timescale, as in |Liu & Showman] ( |2013| l, leads to effectively 
complete cessation of variability as well as sensitivity in three-dimensional (3D) simulations 
with traditional primitive equations. Such momentum (Rayleigh) and thermal (Newtonian) 
drags are, however, ad hoc for 3D giant planet simulations. For 3D hot-Jupiter simulations, 
which typically already employ strong Newtonian drag in the upper region, sensitivity is 
not quenched if only the Newtonian drag is applied in the lower region, without the strong 
Rayleigh drag: in general, both sensitivity and variability persist if the two drags are not 
applied concurrently in the lower region. However, even when the drags are applied concur¬ 
rently, vertically-propagating planetary waves give rise to significant variability in the ~0.05 
to 0.5 MPa region, if the vertical resolution of the lower region is increased (e.g. here with 
1000 layers for the entire domain). New observations on the effects of the physical setup and 
model convergence in ‘deep’ atmosphere simulations are also presented. 
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1 PROLEGOMENON 


Many studies have explored the effects of initial condition and 
resolution or dissipation in hot-Jupiter atmospheric flow simu- 


lations (e.g. Cho et al. 2008| Thrastarson & Cho||2010[ Heng, 

Menou & Phillipps|201 1| Thrastarson & Cho|2011 

|Polichtchouk 

& Cho|2012| Bending, Lewis & Kolb|2012 Liu & S 

iowman|2013 


Polichtchouk et al.||2014| l. |Thrastarson & Cho| j2010| l [hereafter 


TC], in particular, have emphasised the nonlinear effect of initial jet 
configuration on the subsequent evolution, calling attention to it as 
a source of uncertainty for quantitative predictions. ‘Quantitative’ 
here means such things as the precise location of hot regions and the 
three-dimensional (3D) shape and magnitude of vortices and jets, 
which affect the directly-coupled temperature distribution as well 
as the wave momentum and energy depositions. Similar emphasis 
has also been made by |Cho et al.| ( |20(38) l in the two-dimensional 
(2D) context with initial turbulent eddies, in which strength of the 
eddies affected the final mean flow configuration. In this paper we 
revisit the effects with the commonly-used traditional (i.e. hydro¬ 
static) primitive equations with extended, ‘deep’ vertical domain]^ 


* On leave from QMUL; Email: J.Cho@qmul.ac.uk 
^ The discussion here, however, is also germane to understanding results 
from simulations with the compressible, non-hydrostatic, 3D Navier-Stokes 
equations (e.g.|Dobbs-Dixon & Lin|2008[|Mayne et al.|2014|. 


Recently, |Liu & Showman] ( |2013| l [hereafter LS] have per¬ 
formed a set of 3D traditional primitive equations simulations with 
specified thermal (Newtonian) and momentum (Rayleigh) drags 
and assert that TC observe sensitivity because they fail to meet a 
presumptive setup criteria for ‘real’ hot-Jupiters (LS, p. 48)[^In 
their work, LS also advocate the use of a strong Rayleigh drag 
in a fiducial region (denoted as ^ in this paper) at the bottom 
of the modelled atmosphere. Here by ‘strong’ we mean a spec¬ 
ified Newtonian or Rayleigh drag with a timescale of few (or 
less) planetary rotations; and, Si is defined as the pressure range, 
p = [1,20] MP^ where both drags are ‘turned on’ in LS. In 
this paper, we clarify the situation: the insensitivity occurs in LS 
because the combined, strong Newtonian and Rayleigh drags sup¬ 
press effectively all dynamics, everywhere in the domain, except 
for laminar high-speed jets. Quantitative features of the flow, as de¬ 
fined above, are sensitive to initial conditions when the strong drags 
are not applied. 

^ LS also emphasise that the equilibrium temperature and relaxation time 
profiles in TC are constant. However, as stated in TC and explicitly shown 
in |Polichtchouk et S!] (2014j , the sensitivity and variability persist un¬ 
der vertically-varying equilibrium temperature distribution and relaxation- 
timescale that match closely with those used in LS, over the domains con¬ 
sidered in TC and Polichtchouk et al. (2014). 

^ 1 MPa =10 bars 
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Rayleigh drag, which is often used to crudely represent a phys¬ 
ical boundary layer in 3D simulations of large-scale flows (e.g. 
[Held & Suai^[T^ , is implausible for the & region of a giant 
planet without a ‘solid’ surface or with one located deep in its in¬ 
terior - particularly if the drag is a strong one|^The use of such 
a drag in & is notable given that past hot-Jupiter studies without 
the drag have already reported a lack of sensitivity and variability 
(e.g. [Showman et al.||2009[ [Cooper & Showman||2005) : hence, it 
is not required per se to demonstrate insensitivity. Nonetheless, to 
elucidate the effects of the new setup with strong drags, we set up 
here simulations as in LS and cross-check the results with several 
different codes. In doing so, we explicitly demonstrate that the sen¬ 
sitivity is quenched only under the application of an unphysical, 
strong Rayleigh drag and that such a drag is not necessary for at¬ 
taining equilibration if a strong Newtonian drag is already applied[^ 
In any case, ‘equilibration’ reached - with or without the drag - is 
best regarded as heuristic and should not be taken too literally in 
the present setup because unphysical supersonic flow results from 
it. Even with both drags applied, variability and weak sensitivity 
persist if vertical resolution of the ^ region is increased in the sim¬ 
ulations. 


2 SETUP 

TC and LS solve the traditional primitive equations (e.g. [Holton[ 
[2004[ > with the same, or effectively the same, boundary conditions 
(zero ‘vertical velocity’ - i.e. ‘free slip’ - at the top and bottom of 
the domain) [^However, the basic numerical algorithm, treatment 
of explicit viscosity and vertical coordinate are all very different: 
in TC pseudospectral and superviscosity with CAM jCollins et al.[ 
[2004[ > in ? 7 -coordinate and in LS finite-volume and moderate-order 
Shapiro filter with MITgcm ( [Adcroft et al.|20i2) in p-coordinate, 
with the vertical layers equally spaced in log(p). In this study, we 
solve the same equations with the free-slip boundary conditions; 
but, unlike in TC and LS, we employ both the pseudospectral and 
finite-volume models. While the latter is the same model, grid con¬ 
figuration (cubed-sphere) and vertical coordinate used in LS, the 
former is different than the model used in TC: the BOB model 
( [Scott et al.[|2004| in p-coordinates, with vertical layers equally 
spaced in log(p) and p, is used here. BOB is purposely chosen to 
provide an independent check of TC’s results, as well as those of 
LS. Note that, in general, models employing the pseudospectral al¬ 
gorithm converge much faster with resolution than models using 
finite-volume algorithms (see e.g. [Polichtchouk et al.[2014l as well 
as the discussion below). 

In all three studies (TC, LS and the present one), planetary 
parameter values characteristic of the hot-Jupiter HD209458b are 
used (see TC for the values). In this paper, all times and lengths are 

^ Recently, [Schneider & Liu[ (2009} have put forth a weak Rayleigh drag 
representation of ‘magnetohydrodynamics-induced drag’ effect, but the rep¬ 
resentation is at present not widely accepted. 

® The latter point was not made explicit in TC and some confusion appears 
to exist concerning ‘equilibration’ under strong, planetaiy-scale Newtonian 
drag - perhaps unduly influenced by baiotropic non-divergent (i.e. incom¬ 
pressible 2D) turbulence formalism, which does not involve the thermody¬ 
namic fields. 

® MITgcm has an option to allow the bottom pressure surface to be free 
(i.e. nonzero vertical velocity). The state of this option is not specified in 
LS; but, the results are not noticeably affected by the option for the setup 
as in LS. However, in general, energy is added and stability is affected with 
the option on. The option is off in all the simulations discussed in this paper. 


scaled by the planetary rotation period (r = 3.024x10® s) and radius 
{C — 10® m), respectively - unless units are explicitly given. Hence, 
f = 1 corresponds to 1 HD209458b day. Also in all three studies, 
the thermal forcing is crudely represented by a relaxation (drag) to 
a specified ‘equilibrium’ temperature distribution T^[\, (f),p) with 
a drag-time profile tn(p); here A is longitude and </) is latitude. 

Note that and tn are not same in TC and LS, as stressed 
by LS. But, they are not very different: the difference is no more 
than few percent over the model domain in some of the simulations 
reported in TC and all of the forced simulations in [Polichtchouk et| 
[aT] ( [2014) , which uses one of the vertically-varying profiles from 
the TC study. In the present study, the profile Te(p)|A,(#> is also 
not precisely (but essentially) same as that in LS: the profile in LS 
is represented as continuous, piece-wise linear functions in three 
distinct p regions, as tn is in LsQaIso, our Tn in & is 10® s (= 3.3 
scaled), rather than 10^ s, independent of p. The smaller value is 
as in Fig. 6 of LS, and in simulations presented in most of their 
discussion. However, with both drags applied in insensitivity to 
initial flow does not depend on whether the p-independent tn in & 
is 10® s or 10^ s. In fact, as demonstrated below, none of the above 
differences in the physical setup are pertinent to the question of 
sensitivity: only the Rayleigh drag, as employed in LS, is pertinent. 

The crucial Rayleigh drag is characterised by the drag-time 
profile tr{p). Note that, in all the simulations discussed in this 
paper, the Rayleigh drag is applied only in the vertical region, 
log(p) = [5.0, 7.3] = [5.0, 6.0] U^, with tr(p) decreasing expo¬ 
nentially with log(p) - i.e. linearly with p. The specification is as 
in LS. When the drag is not applied in it is not applied anywhere 
in the domain. The strength of the p-dependent tr is characterised 
by min{rR(pG ^)} = rR(20 MPa), and denoted rR 2 o here. 


3 RESULTS 

Fig. □ summarises one of the basic results of this paper: without 
the ad hoc strong drags applied in sensitivity to initial flow 
is robust. The figure shows time-averaged, zonal-mean zonal ve¬ 
locity u* (0, p) from four pairs of simulations with identical setup 
- except for the direction of the initial, barotropic zonal jet (see 
TC for details of the jet profiles) and the presence or absence of 
both Newtonian and Rayleigh drags in In Fig. [T^ and [i]3, "Tn 
and Tr are as in the simulations of LS that apply both drags. In 
Fig. and[^, Tn is as in LS outside but with tn, tr 2 o —> oo 
inside & (i.e. the Newtonian drag is not applied in & and the 
Rayleigh drag is not applied anywhere in the domain). Simulations 
in Fig. [T^ and [TJ; are performed with BOB (T85L40 resolution; 
timestep size. At — 2x 10~'*; and, 8th-order hyperviscosity co¬ 
efficient, izg = 1.5 X 10“^®), and simulations in Fig. and[^ 
are performed with MITgcm using the energy-conserving vorticity 
advection scheme (C64L40 resolution; At = 10“^; Shapiro filter 
order, n = 4; and, filter strength ratio, Al/rghap = 0.25, where Tshap 
is the Shapiro filter coefficient)Here ‘T85L40’ denotes 40 verti¬ 
cal layers with 85 total and 85 sectoral spherical harmonic modes 
in the Legendre expansion per layer, and ‘C64L40’ denotes 40 ver¬ 
tical layers with 6 x 64 x 64 grid points per layer (see [Polichtchouk[ 
[et al.|2014l for details of the models). 

Several features are readily apparent in Fig.[T] First, the East 

^ N.B. Tn here is ‘T^ad’ in LS. 

® Note that MITgcm can be sensitive to the numerical parameter values, as 
we discuss later; hence, simulations that match most closely with the BOB 
simulations are presented. 
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Figure 1. Time-averaged (over t = [250, 300]), zonal-mean zonal velocity u* {<!>, p), where 0 is latitude and p is pressure, from simulations initialised with 
±1000 m s“^ (EastAVest) equatorial jet. Simulations in a) and c) are performed with BOB and in b) and d) are performed with MITgcm in cubed-sphere 
grid. All simulations are performed with vertical layers equally-spaced in log(p [Pa]). In a) and b), tn = 3.3 independent of log(p) and tr in & decreases 
exponentially with log(p) for p > 0.1 MPa with tr2o = 3.3. In c) and d), both tn, tr — ^ oo in 1 ^ (i.e. no drags in &); hence, rR2o oo. When both drags 
are applied, u* is essentially same, in'espective of the initial jet direction [a) and b)]; when both drags are not applied, u* is strongly dependent on the jet 
direction [c) and d)]. 


and West BOB simulations in Fig.^ show no practical difference 
between them, in agreement with LS using the MITgcm. Second, 
not only do the East and West MITgcm simulations in Fig.[T|5 also 
show no practical difference between them, they match very well 
with the corresponding BOB simulations in Fig. [T^. The match is 
not exact, however, and we shall return to this point shortly. For 
now the salient point we wish to make is that the MITgcm results 
here are nearly identical to those presented in hence, there 
is no issue with the aforementioned minor difference in T], or with 
the numerical parameter values chosen for the MITgcm simulations 
in this study. Third, as seen in Fig. [TJ; and [T]i, when the strong 
drags are not imposed in the u* distributions in the East and 
West simulations are clearly distinct - even when the data is heavily 
averaged, as in the figures. Grossly speaking, there are three jets in 
both, East and West simulations, but the width, strength and vertical 
structure of the jets are all noticeably different. Observe as well that 
the East and West simulations in Fig.[TJ; (and Fig. separately do 
not match the corresponding East and West simulations in Fig.[T^ 
(and Fig. [^) - particularly below log(p) « 3 and away from the 
equatorial region. According to the figure, there are at least three 
plausible states, depending on the setup and initial condition. 

The first two features above confirm simultaneously the BOB 
simulations in this study and the MITgcm simulations in LS. The 
very good agreement between the two models (cf. Fig.[T^ an d\^) 
is significant, as this is the first time extrasolar planet atmospheric 
flow simulations from two different numerical models are shown 
explicitly to produce nearly quantitatively same results over the en- 


® cf. with Fig. 8 in LS, modulo plot aspect ratio and minor differences in 
the color range and palette 


tire domain - results which are numerically converged over a good 
range of parameter values. In general, significant variations in the 
simulation results are observed across different models, as well as 
within a single model (see e.g. |Policht chouk e t al.|201^|Thrasta r-| 
|son & Cho|2011| and the discussion below). It is imperative to un¬ 
derstand, however, that the good agreement here is due to the highly 
constraining setup. Nevertheless, as already pointed out, the agree¬ 
ment is still not exact: the equatorial jets in Fig. [TJi are faster and 
sharper than those in Fig.[^, especially outside log(p)« [3.3, 4.3]. 
LS report that in some cases modest variability is exhibited but sen¬ 
sitivity is not exhibited, after the simulation reaches the statistically 
steady state. We have found that this is the case with the setup as in 
Fig.[^ and[^. 

The third feature confirms the behaviour reported in TC: the 
flow is sensitive to the initial state. In particular, it demonstrates 
explicitly that the sensitivity observed in TC does not depend on 
whether the day-night thermal gradient extends all the way down, 
or only part way down, to the bottom of the domain. LS speculate 
that the sensitivity in TC may be caused by a global-scale baro- 
clinic instability induced by the gradient present at the bottom of 
the domain. However, TC have reported sensitivity in a variety of 
setup, including one with reduced equilibrium temperature gradient 
at the bottom. Moreover, |Polichtchouk & Cho| ( |201^ have explic¬ 
itly shown that, if a flow is baroclinically unstable, the instability is 
not thwarted by adding inactive layers (or drags) below the unsta¬ 
ble region. Hence, the placement of the bottom boundary is not the 
crucial factor for sensitivity in these studies. However, such fac¬ 
tors (including numerous other physical and numerical ones) can 
in fact affect instabilities or wave dynamics that naturally arise in 
the course of the flow evolution, steering it to different regions of 
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the configuration (solution) space - and, this is so even if simu¬ 
lations begin with identical initial conditions (e.g.|Polichtchouk &| 
|Cho|201^|Thrastarson & Cho|201 l[ l. The point here is that the state 
illustrated in Fig.[^ (and[Tj5) is an extraordinary one. 

Interestingly, behaviour similar to that in Fig. [T]; and[TJl can 
also be seen in the correspondingly-similar simulations of LS (see 
Fig. 15 therein, bottom rovj^. But, LS interpret the behaviour as 
‘artificial’ because the total angular momentum does not converge 
to the same value. We note here that there is no a priori reason why 
this measure needs to be the same among different model instanti¬ 
ations with different initialisation and balancej^ in addition, even 
if the measure is same, the spatial distribution of the flow (and thus 
temperature) need not be same [^Indeed, this had been precisely 
TC’s point and that of |Cho et al. 1^200^ , as well as the present 
paper: the (unknown) magnitude and distribution of angular mo¬ 
mentum and eddy kinetic energy in the actual atmosphere - needed 
for model initialisation - lead to divergent model evolutions, pre¬ 
venting any kind of precise quantitative predictions at present. 

At this point, we briefly discuss an important numerical issue: 
we have found that MITgcm simulations suffer significant runaway 
angular momentum and acute sensitivity to numerical parameters, 
when both drags are not employed (Fig. Bi)- An extensive study 
of angular momentum conservation in several numerical models 
has been performed and is described elsewhere (Politchouk & Cho, 
in prep.). The basic issue is illustrated in Fig.[^ in which time se¬ 
ries of mass-weighted, global-average, axial angular momentum M 
for simulations with Eastward-jet and Westward-jet initialisations 
at different resolutions are presented^ each series is normalised 
by its initial value, 4.34 x lO^'^kgm^^^ for the East simulations 
and 2.76 x 10®"^ kg m^ s“^ for the West simulations. Fig. [2^ shows 
the runaway and sensitivity (with resolution) in MITgcm simula¬ 
tions. The runaway and sensitivity reduce at high resolution (e.g. 
at C64), but do not vanish. In contrast, runaway and sensitivity are 
not observed in BOB simulations, as seen in Fig. w The runaway 
is detectable in Fig. [T] notice the small north-south hemispherical 
symmetry breaking in Fig. [T]l, but not in Fig. [T];. Since external 
symmetry breaking has not been introduced in these simulations, 
the asymmetry is a manifestation of the runaway - i.e. a numerical 
error. Hence, total angular momentum is not a useful measure for 
MITgcm in this case; and, BOB and CAM results (which are in 
good agreement with each other) cannot be robustly reproduced or 
tested with MITgcm when strong drags are not applied in 

Given the zonal and temporal averaging in Fig. [T] it may be 
difficult to fully appreciate just how different the fields can be in 
simulations initialised with different jets when the strong drags are 
not applied in Unaveraged fields contain more information and 
are more useful for comparing with observations. Fig.j^shows the 
instantaneous flow and temperature fields from the East and West 
simulations of Fig.at the levels, log(p) = {2.0,4.0, 5.6, 7.0}, at 

Note, the Newtonian drag is employed in @ here, unlike in Fig.^; and 

but, the inclusion of the Newtonian drag is not pertinent to the present 
discussion, as shown helow. 

Here ‘balance’ refers to the degree of vortical mode dominance (over 
gravity wave mode) of the flow fie ld (e.g.|Ford, McIntyre & Norton|2000|. 

See |Polichtchouk et al.|(^14} ; this can also be interred from the Fast 
and West pairs of simulations in Fig. [^, presented helow, in which the 
flow and temperature fields are markedly different despite starting with (and 
maintaining) same values of global angular momentum. 

M = (Rp/g) cos (/)-|-tt) cos dU, where O is the rotation 

rate of the planet, Rp is the radius of the planet, g is the gravity, u is the 
zonal velocity, dV is the volume element and ® is the global domain. 



time [t] 



Figure 2. Mass-weighted, global-average, axial angular momentum 
[kgm^ s“^] time series for simulations with East(E)/ West{W) initialisa¬ 
tions at different resolutions. Each series is normalised to its initial value: 
4.34 X lO^^kgm^ s“^ (E) and 2.76 X 10^‘^kgm^ s“^ (W) simulations. 
The physical setup is as in Fig.[Tjl (MITgcm) and Fig.[T}; (BOB). MITgcm 
exhibits monotonic runaway behaviour, which reduces with increased res¬ 
olution. In contrast, the BOB model conserves angular momentum exactly, 
independent of the resolution. 

t — 300. Both the flow and temperature distributions are markedly 
different in the two columns (East and West simulations), except 
near the top of the domain (cf. log(p) = 2 frames). The behaviour 
near the top is expected, given the extremely strong Newtonian drag 
applied there (tn « 0.03). But, even there significant differences 
are present in both fields, as will be shown more clearly below. 

Significantly, there is much emphasis in the current literature 
on the temperature distribution at high altitude (e.g. p = 30 mbar 
=> log(p) ~ 3.5), where both the temperature and the flow are 
highly restrained by the very short Tn (« 0.1). However, attention 
should also be given to the lower regons - for a more realistic, com¬ 
plete picture!^ First, note that such timescales are comparable to 
the Brunt-Vaisala period, the timescale of ‘deep’ internal gravity 
waves. The distortion or omission of various waves, including the 
gravity waves, is a source of significant inaccuracy in current at¬ 
mospheric flow simulations: for example, wave interactions with 
the background flow (which modify the temperature distribution) 
are poorly captured - if at all (see e.g. |Ford, McIntyre & Norton) 

This also applies to the mid- and high-latitude regions. 
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Figure 3. Instantaneous latitude-longitude maps of the flow v (m s~^) and temperature T (K) fields at four pressure levels, log(p) = {2.0,4.0, 5.6, 7.0}, at 
t = 300 in cylindrical-equidistant projection from the {West/East) simulations in Fig.^; v = {u, v), where u and v are the zonal and meridional velocities, 
respectively. The reference flow vectors are shown at the bottom right in each panel and temperature ranges for each row are shown at the right. In both West 
and East simulations, the flow and temperature distributions are complex with multiple, irregular hot/cold regions - often situated away from the equator 
and substellar/antistellar (0°/180°) longitude. The spatial distribution of the flow and temperature fields are significantly different, except near the top of the 
modelled domain (see Fig.[^. Symmetry has not been broken initially, unlike in most simulations in TC. The vertical temperature gradient fields in the two 
simulations are also different, and would lead to different emergent heat flux distributions. 


|200Q[ |Cho et al.|2Q03[ [Watkins & Cho||20I0[ and also this study 
below). These interactions can be long-range and be effected by 
waves originally of small amplitude or scale. Second, infrared flux 
itself can originate from any of the levels modelled, consistent with 
the use of Newtonian relaxation approximation. The approximation 
is meant to represent crudely the effect of radiation emanating from 


the atmospheric region where the approximation is being invoked 
(e.g. |Salby|p^96[ [Andrews et al.|p~987[ [Cho et ^|2008| t. Given 
the above, matches with observations in current simulations may 
merely be fortuitous: a better understanding of both the physics 
and numerics is needed. 

Note also that the fields in Fig. are very similar to those in 
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Figure 4. West — East difference maps of the zonal wind u (m s“^) and temperature T (K) fields at the p-levels in Fig. The wind differences are at the 
left column, and the temperature differences are at the right column. The plot ranges ai'e shown at the right of each panel. Locally, the absolute flow and 
temperature differences can be large - as much as ~4000 m s“^ and ~200 K at log(p) = 2, respectively. 


TC, despite the bottom boundary being located at a much greater 
depth in the simulations of the figure. As in TC, the hottest and 
coldest regions are not simply connected (e.g. |Mutrkres|2000l l and 
often located away from the equator. This is in contrast to ‘maps’ 
constructed from observations which assume a latitudinal distribu¬ 
tion monotonically decaying away from the equator. Moreover, the 
hottest (coldest) region can often be situated at the night (day) side, 
also in contrast to what has been reported in some observations and 
simulations (e.g. |Knutson et al.|2007| LS, and references therein). 
In general, the locations of the temperature extremum regions vary 
in space and time, within a single simulation and across different 
simulations (e.g. |Cho et al.|200^|2008[ [Thrastarson & Cho|2010[ 


IPolichtchouk et al.|2014) , again depending on the initial flow and 
equilibrium temperature states specified. 

For a clearer picture of the magnitude and spatial distribution 
of the dijferences. Fig. presents point-wise subtractions of the 
West and East frames at the p-levels shown in Fig.|^ In Fig.|^ the 
left column shows the instantaneous zonal wind field difference, 
Uw(A, (f)) — Ue(A, (j)), and the right column shows the instantaneous 
temperature field difference, Tfi{X, 4>) — 2e(A, (f>)- Note the large, 
absolute maximum flow and temperature differences (more than 
~4000 m s“^ and ~200 K, respectively) locally at log(p) = 2. 
At greater depths the u-differences are reduced, but are still quite 
large (~1000 m s“^) and on global scales. The u-difference in 
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the left bottom panel reflects the opposite sign of the zonal jets in 
the ^ region in the two simulations, particularly in the equatorial 
region (cf. East and West simulations in Fig.[TJ;, for example). The 
T-differences also generally decrease with depth. Notably, in the 
3) region, the temperature difference is a very small fraction of the 
ambient temperature (< 4%) in both of the simulations, indicating 
that the Newtonian drag as specified in simulations of Fig.[^ and 
[T]} is not really effective or needed. This is not surprising, given the 
large inertia of 

In their study, LS motivate the use of a strong Rayleigh drag 
in as an ‘accelerator’ to help reach the ‘equilibrium’ state more 
quickly. However, such a drag can steer the flow to artificial states 
and is in actuality unnecessary - particularly if temperature is used 
to characterise the equilibration, as in many climate studies (e.g. 
[Brandefelt & Otto-Bliesn^|2009[ >, or if strong Newtonian drag 
is already applied in & (which is also dubious, as noted above). 
Both can be clearly seen in Fig. where time series normalised 
by the initial value from four long-duration simulations (roughly 
‘permutations of {( 71 ( 20 ,Tn)}’) are presented; all the simulations 
here employ strong Newtonian drag in the upper region. In the fig¬ 
ure, the mass-weighted, global-average temperature (T) is steady 
in all the simulations presented, regardless of the drags employed 
in tSi (black curve). With only the Newtonian drag applied in the 
mass-weighted, global-average kinetic energy {X) reaches steady 
state (at t « 700, red curve). With both drags not applied in 
slow increase in {X) is observed (green curve), but the rate of in¬ 
crease depends on time (and on the model setup). We stress that, in 
principle, such an evolution is physically valid - as long as {tX) 
does not ‘blow up’. Observe that the Rayleigh drag is responsible 
for suppressing the sensitivity (cf. orange, blue and red curves): for 
this, Newtonian drag by itself is ineffectual. 

FS also suggest, following e.g. the study by|Perna, Menou| 
|& Rauscher| ( |2010l l for higher altitude, that a Rayleigh drag in @ 
might serve as a crude representation of ‘magnetic drag effects’ 
stemming from thermal ionisation. There are two major concerns 
with this. First, thermal ionisation is insignificant in the ^ region: 
temperature is too low and density is too high. This is so even tak¬ 
ing into account the low ionisation potential of alkali metals (e.g. 
K, Na, Ca), as these are trace species (see e.g. |Lewis|2004^ . Using 
solar abundances, njj+ /rin < 2 x 10“^® andn^r- /njj+ ~ 3x10®, 
where is the a:-specie number density and ‘n’ subscript refers to 
neutral. Hence, the bulk ionisation level (electron volume mixing 
ratio), Xe = rLelrin^ ~ nj,^+ is at least 10 ^ times lower than 
that required for the fluid medium to be influenced electromagnet- 
ically: Xe ~ 10 “^ is required for ionisation drag effects to become 
significant (e.g. |Schunk & Nagy|2000[[Koskinen et al.|2014^ . Sec¬ 
ond, even if the ion-induced drag were significant via a non-thermal 
mechanism, it cannot be represented as an isotropic drag to rest 
on the momenturrF^ ion velocities and the intrinsic field orienta¬ 
tion need to be modeled self-consistently for accurate representa¬ 
tion (e.g. |Koskinen et al.|2010^ . 

Fet us now consider a more detailed look at the behaviour 
in time. Fig. [^presents t-p Hovmdller plots of the zonal velocity, 
w(A = 0° 0 = 30°, p, f),forf = [200, 300] and log (p) = [4.0, 7.3]. 
Fig. 1^ and 1^ are from the simulations in Fig. which is at 
T85L40 resolution. Their general behaviours have been verified 
with up to T341L40 resolution using the BOB model. Fig. 
presents a simulation with a setup identical to that of the simu¬ 
lation in Fig. [T^ (East), except the resolution here is T21F1000 


This is also true of gravity wave induced drag. 
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Figure 5. Mass-weighted, global-average temperature (T) (black) and ki¬ 
netic energy {X) (green, red, blue and orange) time series from four sim¬ 
ulations set up identically, except for different (tr 2 o, Tn) in St. When finite, 
tr(p) decreases linearly with p in log(p) = [5, 7.3] = [5, 6] and tn 
is constant in S\ outside their respective regions, rN(p) and tr(p) are as 
in LS in all cases. All are T85L40 BOB simulations, initialised with the 
-1-1000 m s“^ jet and To = 1800 K. The black curve is common to all four 
simulations. The blue and orange curves are very close to each other. The 
green and red curves are different for a different initialisation while the blue 
and orange curves are not. 


(Af = 8 X 10“"^ and tzg = 3 x 10“^®). The results are nearly 
identical at T42L500 and T85L40 resolutions. The simulation at 
T21L1000 resolution is presented here for equatable comparison 
with the simulation shown in Fig. 1^, in which the setup is identi¬ 
cal to the simulation shown in Fig.lM, except the vertical layers are 
now equally spaced in p. This spacing has the effect of representing 
the lower (higher) region with greater (fewer) number of layers: the 
resulting resolution in Fig. for example, is 57 times greater at 
the bottom of ® than that in Fig.|^ and in LS. In the simulation 
presented in Fig.|^, the spacing is more numerically consisten^^ 
as the equations solved in all the simulations are in p rather than 
log(p) coordinate. 

Fig.|^ and|^ show vertically-propagating Rossby (planetary) 
waves, generally prominent in the model atmospheres when strong 
drags are not applied in &. In the two panels, the waves are propa¬ 
gating in opposite directions (upward and downward in Fig. |§i and 
respectively), as evident from the tilt of the phase lines. The 
magnitude of the peak amplitude is ~700 m s“^ in both, given the 
latitude location of the constructed plots, but the sign and period 
are different - as is the growth with height over time. As already 
discussed, such waves can induce significant modification of the 
background flow via saturation and encounter with critical layers 
(e.g. [Andrews et al.|198'7l |Holton|2004[ ). The action of the waves 
will be discussed more in detail elsewhere. Here we wish to high¬ 
light the suppression of essentially all temporal activity - including 
these waves - when the strong drags are applied in in addition 
to a strong Newtonian drag already applied throughout nearly all of 
the domain outside ^ (cf. Fig. 1 ^ with Fig.|^). While not uninter¬ 
esting, such a ‘dead’ atmosphere shown in Fig. [5 is not dynamic, 
and a flow simulation of it is not very informative. 

Interestingly, the suppressed behaviour is not robust and ap¬ 
pears to be a numerical artifact. This is demonstrated in Fig. B- 


Numerical consistency refers to discretisation error tending to zero as 
the resolution is increased, under pointwise convergence at each grid point 


(e.g. |Strikwerda|2Q0^ . 
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Figure 6. t-p Hovmoller plots of instantaneous zonal velocity 3 q) for t = [200, 300] and log(p) = [4.0, 7.3]. The colour bars indicate 

the velocity for each row. All simulations are performed with BOB. Initial jet amplitudes and directions (indicated above each plot) are as in Fig.^ In a) 
and b), the resolution is T85L40 and tr 2 o, tn ^ oo in In c) and d), the resolution is T21L1000 and rR 2 o, tn 3.3 in ^ (as in LS) with the vertical levels 
equally-spaced in log(p) (c) and p (d). Vertically-propagating waves of planetary scale power variability and behave differently, depending on the initial flow, 
when strong drags are not applied in ^ [cf. a) and b)]. However, even with the strong drags employed, propagating waves are present if the vertical resolution 
is increased in ^ [cf. c) and d)]. 


The panel shows the behaviour of the simulation in Fig. when 
the vertical resolution of the dynamically active region (log(p) ~ 
[5, 6]) is increased, principally by employing a different spacing. 
As can be seen in Fig. |^, the variability returns - powered by 
waves that strongly shear, emerging from the & region, and in¬ 
duce secondary and even tertiary waves upon interacting strongly 
with the northern flank of the equatorial jet at log(p) ~ [5.0,4.5] 
(see Fig.[T^). The amplitude of the waves here is much smaller than 
those in Fig. E an d[6j5, consistent with the much stronger drag in 
this case. Propagating gravity waves that are similarly generated 
can enhance the observed variability as well (e.g. |Scinocca & Ford| 
|2000[ [Watkins & Cho|2010| >, but these waves are poorly captured 
in this simulation at the employed resolution. 


4 CODA 

We close with some thoughts and observations arising from this 
study. We have carefully cross-checked the results with codes used 
in TC and LS, as well as a third code (which has been extensively 
tested under conditions appropriate for hot-Jupiters). As shown by 
[Polichtchouk et aL| ( [2014^ and the present study, one can arrive at 
erroneous conclusions if simulations from different codes are not 
at least qualitatively reproduced with the same setup. Even then, 
erroneous conclusions can still be drawn, as codes which perform 
well ostensibly in one region of the physical parameter space do not 


perform well in another (or, more precisely, extended) region. For 
example, when pushed to a highly ageostrophitf^region (as would 
occur in a typical hot-Jupiter atmosphere simulation), numerical 
accuracy of the code can become seriously degraded by small-scale 
oscillations generated in that region of the parameter space (e.g. 
[Thrastarson & Cho|201 1^ . 

By expanding the study by [Polichtchouk et al.[j2014) into the 
initial condition space, the results here constitute an extension of 
that study: in general, simulations are sensitive to initial condition - 
in addition to many other factors, physical and numerical. That the 
physical system under study exhibits multiple equilibrium states 
(hence, also non-ergodicity) is not, it seems to us, a particularly 
startling assertion given its highly nonlinear and poorly constrained 
forced-dissipative nature. But, it has been challenged forcefully by 
LS. 

Ultimately, it appears the question of sensitivity and variabil¬ 
ity rests on whether one believes that it is truly ‘realistic’ to apply 
Rayleigh drag as well as Newtonian drag on very short timescales, 
particularly in the & region. Even assuming that the drags are 
acceptable representations of the forcing and dissipation in hot- 
Jupiter atmospheres, one must ask; how realistic are tn, Tr and Te 
currently used? The responses of TC and LS clearly differ on this 
question. A telling observation of this study is the extraordinary 


e.g. high Rossby and Fronde numbers - measures of rotation rate and 
stratification, respectively (see e.g. [Holton[2004) 
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length one must go to suppress sensitivity and variability inherent 
in the system: all eddies and waves must be eliminated for all time. 
One may wonder whether such a state is really realistic, given that 
the atmosphere would likely contain some vorticity and turbulence 
(both inhomogeneous and intermittent) and that there are altitude 
regions (e.g. just above characterised by neither short Newto¬ 
nian nor short Rayleigh drag timescales. In our view, the general 
question of realistic forcing and initialisation - and their effects - 
is still unsettled, and the question deserves much more attention 
and scrutiny than has generally received thus far. 

Most importantly, the investigation should proceed with the 
following observations squarely in the fore. We have found that 
MITgcm in the cubed-sphere grid configuration conserves angular 
momentum poorly and is acutely sensitive to numerical parameter 
values in simulations initialised with high speed Jets, particularly 
without the strong Rayleigh drag applied in & (e.g. Fig.|^: hence, 
it cannot be used for assessing sensitivity in this case, as it is not 
possible to establish a reliable baseline for quantification. We have 
also found that, all codes used in this study (including BOB) lead 
to a zonally-symmetric, superrotating, supersonic jet at the equa¬ 
tor, at log(p) « 3. Given that the hydrostatic primitive equations 
with free-slip boundary conditions (as well as with the viscosity 
representation for small Mach number flows) are solved in all the 
numerical models in this study, such flow is physically not valid 
[see e.g. discussion in |Holton]p004| >]. Therefore, claims of realism 
in these setups is moot: such jets may occur on real hot-Jupiters, 
but not in these simulations. Either a different setup must be used 
for the equations and boundary conditions solved or a different set 
of equations and boundary condition must be solved for the setup 
used. Relatedly, while all of the simulations presented here exhibit 
supersonic zonally-symmetric equatorial flow at some height, not 
all simulation do: in general, that depends on the physical setup, 
numerical scheme and initial condition employed (e.g.|Thrastarson| 
|& Cho|2010||Polichtchouk et al.|2014| l. 


ACKNOWLEDGMENTS 

The authors acknowledge helpful discussions with Craig Agnor, 
Tommi Koskinen and Stephen Thomson - and particularly for their 
careful reading of the manuscript. The authors also acknowledge 
the hospitality of the Kavli Institute for Theoretical Physics, Santa 
Barbara, where some of this work was completed. I.P. is supported 
by UK’s Science and Technology Facilities Council research stu¬ 
dentship. H.Th.Th. is supported by the NASA Post-doctoral Pro¬ 
gram at the Jet Propulsion Laboratory, administered by Oak Ridge 
Associated Universities through a contract with NASA. 


Cooper C. S., Showman A. R, 2005, ApJ, 629, L45 
Dobbs-Dixon, I., Lin, D. N. C., 2008, ApJ, 673, 513 
Ford, R., McIntyre, M. E., Norton, W., A. 2000, J. Atmos. Sci., 
57, 1236 

Held I. M., Suarez M. J., 1994, Bull. Amer. Met. Soc., 75, 1825 
Heng, K., Menou, K., Phillipps, P, 2011, MNRAS, 413, 2380 
Holton J. R., 2004, Introduction to Dynamic Meteorology, 4th ed.. 
Academic Press, San Diego 

Koskinen T, Cho J. Y-K., Achillios N., Aylward A. D. 2010, ApJ, 
722, 178 

Koskinen T, Yelle, R., Lavvas, P, Cho J. Y-K. 2014, ApJ, 796, 16 
Knutson H. A., Charbonneau, D., Allen, L. E., Fortney, J. J., Agol, 
E., Cowan, N. B., Showman, A. R, Cooper, C. S., Megeath, S. 
T., 2007, Nature, 447, 183 

Lewis J. S., 2004, Physics and Chemistry of the Solar System, 2nd 
ed.. Academic Press, Amsterdam 
Liu B., Showman A. P, 2013, ApJ, 770, 42 
Mayne,N. J., Baraffe, I., Acreman, D. M., Smith, C. Browning, 
M. K., Amundsen, D. S., Wood, N., Thubum, J., Jackson, D. R., 
2014, A& A, 561,24 

Munkres J. R., 2000, Topology, 2nd ed., Prentice Hall, Upper Sad¬ 
dle 

Perna, R., Menou, K., Rauscher, E., 2010, ApJ, 719, 1421 
Polichtchouk, I., Cho Y-K. J., 2012, MNRAS, 424, 1307 
Polichtchouk, I., Cho Y-K. J., Watkins C., Thrastarson H. Th., 
Umurhan O., Juarez, M, T, 2014, Icarus, 229, 355 
Salby M. L., 2000, Eundamentals of Atmospheric Physics, Aca¬ 
demic Press, San Diego 
Schneider T, Liu J., 2009, JAS, 66, 579 

Schunk, R. W., Nagy, A. E. 2000, Ionospheres: Physics, Plasma 
Physics, and Chemistry, Cambridge University Press, Cam¬ 
bridge 

Scott R. K., Rivier L., Loft R., Polvani L. M, 2004, NCAR Tech¬ 
nical Note No. 456 

Scinocca J. E., Eord, R. 2000, J. Atmos. Sci., 57, 653 
Showman A. R, Fortney J. J., Lian Y, Marley M. S., Freedman R. 
S., Knutson H. A., Charbonneau D., 2009, ApJ, 699, 564 
Strikwerda, J. C., 2004, Finite Difference Schemes and Partial 
Differential Equations, 2nd ed., SIAM, Philadelphia 
Thrastarson H. Th., Cho J. Y-K., 2010, ApJ, 716, 144 
Thrastarson H. Th., Cho J. Y-K., 2011, ApJ, 729, 117 
Watkins C., Cho J. Y-K., 2010, ApJ, 714, 904 


REFERENCES 

Adcroft A. et al., 2012, ‘MITgcm User Manual’ 

Andrews, D. J., Holton J. R., Leovy, C., 1987, Middle Atmosphere 
Dynamics, Academic Press, San Diego) 

Bending V. L., Lewis S. R., Kolb U., 2012, MNRAS, 428, 2874 
Brandefelt J., Otto-Bliesner B. L., 2009, GRL, 36, L19712 
Cho J. Y-K., Menou K., Hansen B. M. S., Seager S., 2003, ApJ, 
587, LI 17 

Cho J. Y-K., Menou K., Hansen B. M. S., Seager S., 2008, ApJ, 
675,817 

Collins W. D. et al., 2004, Description of the NCAR Community 
Atmosphere Model (CAM 3.0) 


© yyyy Ras, mnras ooo,[T]-?? 








